Collisionless relaxation in non-neutral plasmas 
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A theoretical framework is presented which allows to quantitatively predict the final stationary 
state achieved by a non-neutral plasma during a process of collisionless relaxation. As a specific 
application, the theory is used to study relaxation of charged-particles beams. It is shown that a 
fully matched beam relaxes to the Lynden-Bell distribution. However, when a mismatch is present 
and the beam oscillates, parametric resonances lead to a core-halo phase separation. The approach 
developed accounts for both the density and the velocity distributions in the final stationary state. 



Relaxation to a final stationary state (SS) of parti- 
cles interacting through long-range forces, such as (un- 
screened) Coulomb or gravitational, is intrinsically differ- 
ent than that of systems with short-range interactions [l[ • 
In the latter case, the interparticle collisions drive the 
system to an equilibrium state described by the Maxwell- 
Boltzmann distribution. This distribution is unique, in 
the sense that it is completely determined by the glob- 
ally conserved quantities such as the total energy, mo- 
mentum, angular momentum, etc. — and is otherwise 
independent of specific initial conditions. The same is 
true for neutral plasmas for which the bare Coulomb po- 
tential is dynamically screened by the counter-charges, 
leading to a well defined thermodynamic limit and equi- 
librium 2]. Relaxation of particles interacting by long- 
range (unscreened) potentials, on the other hand, is very 
different. For these systems, the collision duration time 
diverges and the state of thermodynamic equilibrium is 
never reached. Instead, the dynamics evolves to a sta- 
tionary state in which distribution functions appear to 
stop varying with time. Unlike thermodynamic equilib- 
rium, in SS, however, detailed balance is violated [3j and 
neither equilibrium thermodynamics nor equilibrium sta- 
tistical mechanics can be used. 

In the limit in which the number of particles goes to 
infinity (N — > oo) while the total mass and charge are 
fixed, the non-neutral plasma is described exactly by the 
Vlasov equation 0], 
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where D is the advective derivative, /(t, r, v) is the one 
particle distribution function, and F is the mean force 
felt by particles at position r. For simplicity we have 
set particle mass to unity. Vlasov equation shows that 
the distribution function evolves in time as an incom- 
pressible fluid. If we now discretize the height of the 
initial distribution function / (r,v) into a set of lev- 
els 7/j, with j = l...p, the Vlasov dynamics of a d 
dimensional system preserves each level's hypervolume 
7(^0 = / ^(/(^ r i v ) — ?7j)d d rd d v. For a general dis- 
tribution function this condition is equivalent to exis- 
tence of an infinite number of dynamical invariants called 
the Casimir integrals or simply the Casimirs [5j]. One of 
the Casimirs is the Boltzmann entropy which is, there- 
fore, a constant of motion. While the fine-grained dis- 



tribution function f(t,r,v) never reaches a stationary 
state - the evolution continues on smaller and smaller 
length scales ad infinitum - Lynden-Bell argued that 
the coarse-grained distribution function /(£, r, v, ), aver- 
aged on microscopic length scales, will rapidly relax to a 
meta-equilibrium with /(r, v). For gravitational systems 
Lynden-Bell called this process "a violent relaxation" [f| . 
To obtain the stationary distribution /(r,v), we divide 
the phase space into macrocells of volume d d r d d v, which 
are in turn subdivided into v microcells, each of volume 
h d . As a consequence of incompressibility, each microcell 
can contain at most one discretized level r/j . The number 
density of the level j inside a macrocell at (r, v) — num- 
ber of microcells occupied by the level j divided by v — 
will be denoted by pj(r,v). Note that by construction, 
the total number density of all levels in a macrocell is 
restricted to be 
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r,v) < 1 
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Using a standard combinatorial procedure Q it is then 
possible to associate a coarse-grained entropy with the 
distribution of {pj}. The entropy is found to be that of 
a p species lattice gas, 



S = -k B J ^ij^^.^^ln^^v)] (3) 

+ [l-^ Pi (r,v)]Ml-^ ft (r,v)]l , 
j=i 0=1 J 

where ks is the Boltzmann constant. Lynden-Bell ar- 
gued that collisionless relaxation should lead to the den- 
sity distribution of levels which is the most likely, i.e. 
the one that maximizes the coarse-grained entropy, con- 
sistent with the conservation of all the dynamical in- 
variants — energy, momentum, angular momentum and 
the hypervolumes j(r]j). In terms of the number den- 
sities {pj} which maximize the coarse-grained entropy, 
the stationary distribution function becomes /(r, v) = 
T)jPj(T,v). The maximum entropy state, however, 
can only be achieved if there is a sufficient ergodicity 
(mixing) in the phase space. 
If the initial distribution fo(r,v) is uniform (p = 1), 
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the maximization procedure is particularly simple, yield- 
ing a Fermi-Dirac distribution, 

/(r, v) = ffc P (r, v) = e/3[e(r<v) _ M] ~ 1 , (4) 

where e is the mean energy of particles with velocity v at 
position r, and (3 and \x are the two Lagrange multipliers 
required by the conservations of energy and number of 
particles, 

|d d rd d v £ (r,v)/(r,v)= eo (5) 

J d d rd d v/(r,v) = 1 . 

In the above formula eo is the energy per particle spec- 
ified by the original distribution fo- For an azimuthally 
symmetric system, the mean particle energy e is a func- 
tion of only the modulus r and v. By analogy with the 
usual Fermi-Dirac statistics, we define (3 — 1/fcsT, were 
T is the effective temperature of the stationary state (not 
to be confused with the usual definition of temperature 
in terms of the average kinetic energy which is valid only 
for classical systems in thermodynamic equilibrium) and 
/x is the effective plasma chemical potential. 

In this Letter, we will show that when applied to 
non-oscillating confined non-neutral plasmas, Eq. ([4]) de- 
scribes very accurately the final stationary state. On 
the other hand, if during the relaxation dynamics plasma 
undergoes collective oscillation, the theory of violent re- 
laxation fails dramatically. Instead, we observe that the 
system separates into two coexisting phases — a cold core 
surrounded by a halo of highly energetic particles. The 
relaxation process is extremely slow, taking tens of thou- 
sands of plasma oscillations to reach the stationary state. 
A new approach will then be presented which quantita- 
tively predicts the phase-space distribution functions in 
the final relaxed state. 

To illustrate the general theory, we will apply it to 
study the transport of intense, continuous, charged- 
particles beams through a uniform focusing magnetic 
field [8| . The beam is assumed to propagate with a con- 
stant axial velocity v z e z , so that the axial coordinate 
s — z — v z t plays the role of time. The external focusing 
field is given by B = B Q e z and is used to compensate the 
repulsive Coulomb force between the beam particles. It is 
convenient to work in the Larmor frame of reference Q , 
which rotates with respect to the laboratory frame with 
angular velocity Q,l = qB /2~fbmc, where c is the speed 
of light in vacuo, and q, to, and jb = [1 — (v z / c) 2 ] -1 / 2 
are the charge, mass, and relativistic factor of the beam 
particles, respectively. In this frame, the beam distribu- 
tion function fb(s, r, v) evolves according to the Vlasov- 
Poisson system [8] 

|j + v • Vfb + (-k„t - VVO • V v / fc = 0, (6) 
V 2 = -(27tK/N b )n b (r,s), (7) 



where N b is the number of particles per unit axial 
length, r is the position vector in the transverse plane, 
and v = dr/ds is the transverse velocity, n&(r, s) = 
Nb J fb d 2 v is the transverse beam density profile, k z = 
q 2 B 2 j\^ 2 v 2 z m 2 c 2 is the focusing field parameter, and 
K = 2q 2 Nb I '^v 2 mc 2 is the beam perveance, which is 
a measure of the beam intensity. In Eqs. ^ and 0, 
-0 is a scalar potential that incorporates both self-electric 
and self-magnetic fields, E s and B s . We shall take zero of 
the scalar potential to be at r w , the position of the con- 
ducting channel wall. The distribution function is nor- 
malized, so that J/f,d 2 rd 2 v = 1. In the Larmor frame, 
the system corresponds to a two dimensional non-neutral 
plasma of pseudo particles of mass m p = 1 and charge 
q = y/ K/Nb interacting by a repulsive logarithmic po- 
tential ip(r) = —q 2 hi(r /r w ), confined in a parabolic po- 
tential well of U(r) = K z r 2 /2. We will now explore the 
relaxation of these particles from the initially uniform 
distribution (p = 1), 

/o(r, v) = m Q(r m - r)e(v m - v) (8) 

with rji — l/7r 2 r 2 re i^ l , to the final stationary state. 

At time t = the particles are uniformly distributed 
in the phase space between r < r m and p < p m . The 
distribution function Eq. ([5]), however, is not a station- 
ary solution of the Vlasov-Poisson system, and for t > 
the system will start evolve in time. It is possible to 
adjust the values of r m and v m in such a way that dur- 
ing the evolution, the beam envelope (rms particle posi- 
tion) oscillates as little as possible. This corresponds to 
the so called matched beam condition — the beam re- 
laxes to equilibrium, but without undergoing significant 
macroscopic oscillations. For the distribution function 
([8]) , the matching condition can be determined using the 
beam envelope equation @, Q. It is possible to show 
that the beam will oscillate only little if v 2 n w c z r^ — K. 
When this condition is met, we expect the mixing to be 
efficient and Lynden-Bell theory to apply. The coarse- 
grained beam distribution should then relax to Eq. ([][]), 
with e(r, v) = v 2 /2 + U(r) + tp(r), where the mean elec- 
trostatic potential ip(r) is determined self-consistently by 
an iterative solution of Eq.([7]), subject to constraints of 
Eqs. iJSJl with energy per particle given by 

To compare the theory with the simulations, we calcu- 
late the number particles inside shells located between 
r and r + dr, N(r)dr = 2irNbrdr J d 2 v/(r, v); and the 
number of particles with velocities between v and v + dv, 
N(v)dv — 2nNbvdv J d 2 r/(r, v). In Fig. 1 the solid lines 
show the values of N(r)/Nb and N(v)/Nb obtained using 
the theory described above, while points are the result of 
a self-consistent N-particle dynamics simulation [9(. In 
all the figures distances are measured in units of \J Kj k z 
and velocities in units of \f~K. Excellent agreement be- 
tween the theory and the simulation is found for both 
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position and velocity distributions without any fitting pa- 
rameters. We have checked that agreement persists for 
other values of r m and v m , as long as the matching condi- 
tion is satisfied. The agreement, however, disappears as 
soon as the matching condition is violated and the beam 
begins to oscillate, Figs. (|2|3|) . 




FIG. 1: Position and velocity distributions for a matched 
beam with r m — 1.48a/A7k z and v m = l.l^K. Solid line is 



i^J K/k z and v m = l.lVK. 
the theoretical prediction obtained using distribution function 
of Eq. (j4| , while points are the result of dynamics simulation 
with N = 5000 particles. 

Plasma oscillations lead to a number of important con- 
sequences which are not taken into account in the the- 
ory of violent relaxation. For space-charge dominated 
inhomogeneous beams, the oscillations result in propa- 
gating density waves which eventually break, emitting 
high energy particle jetsUJj] . The oscillations also excite 
parametric resonances transferring large amounts of 
energy to some particles at the expense of the rest |9j, 
see Fig. 2. Both of these mechanisms lead to inefficient 
phase space mixing and non-ergodicity. As the relax- 
ation proceeds, the oscillating beam core becomes pro- 
gressively colder, while a halo of highly energetic parti- 
cles is created. Because of the incompressibility restric- 
tion imposed by the Vlasov dynamics Eq. ((2]), the core, 
however, can not freeze - collapse to the minimum of 
the potential energy. Instead, the distribution function 
of the core particles progressively approaches that of a 
fully degenerate Fermi gas. 

The extent of the halo is determined by the location of 
the parametric resonance, and its range can be calculated 



using the canonical perturbation theory . In Fig. 2a 
we show the Poincare section of a test particle moving un- 
der the action of an oscillating beam potential calculated 
using the envelope equation [{J, [HJ • The resonant orbit is 
the outermost curve of the Poincare plot. The first reso- 
nant particles move in almost a simple harmonic motion 
with energy e R = —Khx{r R /r w ) + K z r R /2, where r R is 
the intersection of the resonant trajectory with the v — 
axis. As more and more particles are ejected from the 
beam core their motion, however, becomes chaotic and 
a halo distribution becomes smeared out. We find that 
the distribution function of a completely relaxed halo is 
very well approximated by the Heaviside step function 
e(e R -e). 

For an out of (thermodynamic) equilibrium system, 
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FIG. 2: (a) Poincare section of a test particle moving in 
an oscillating potential controlled by the envelope equation. 
The outermost curve shows the resonance which determines 
the halo's range, (b) A snapshot of the phase space par- 
ticle distribution obtained using dynamics simulations with 
N = 5000, after 40 thousand beam envelope oscillations. The 
halo particles are almost uniformly distributed in the phase 
space. The halo's range vr is determined by the original para- 
metric resonance, see panel (a), (c) Position and (d) velocity 
distributions. Solid curves are the theoretical predictions ob- 
tained using the distribution function of Eq. (|10[) and points 
are the results of dynamics simulations. The initial distribu- 
tion is uniform with r m = 1.98-y/ K/k z and v m = 0.24VA 7 . 
It has exactly the same energy as the fully matched distribu- 
tion of Fig[TJ showing that SS depends explicitly on the initial 
distribution. 



there are no clear parameters which will control the core- 
halo coexistence — such as pressure, temperature and 
chemical potentials for usual thermodynamic systems in 
coexistence. We can not, therefore, a priory say when 
the halo formation will stop and a stationary state be es- 
tablished. Empirically, however, we have observed that 
this happens when the core temperature becomes suffi- 
ciently low. In all cases studied, we find that the core- 
halo equilibrium is achieved when the ratio between the 
core temperature and the corresponding Fermi tempera- 
ture is T/T F » 1/40 - i.e. when /3/x w 40. The distribu- 
tion function for the core-halo system, then, takes a very 
simple form 



A(r,v) 



3 /3e(r,v)-40 i 1 



(10) 



Since all the dependence on r and v enters only implicitly 
through e, /& automatically satisfies the Vlasov-Poisson 
system. The value of rji — l/^r^v^, is determined by 
the initial distribution /o, while the value of e R is cal- 
culated from the location of the parametric resonance, 
Fig. 2a. This leaves to determine self-consistently, using 
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Eqs. m and (JT]), the mean electrostatic potential ifj(r), 
the inverse temperature f3, and the amplitude \ which 
will determine the fraction of particles inside the halo. 
These can, once again, be obtained iteratively. In Figs. 
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FIG. 3: (a) Position and (b) velocity distributions. Points are 
the result of dynamics simulations. Solid curves are the the- 
oretical predictions obtained using the distribution function 
of Eq. (|10[). Dashed curves are the predictions of the violent 
relaxation theory based on Eq. Q. The figure demonstrates 
that for oscillating beams, mixing is inefficient and violent 
relaxation theory does not apply. The initial distribution is 
uniform with r m = l.Ov/ K/k z and v m = 2Ay r K. 

[2] and [3] we plot N(r)/N b and N(v)/N b , obtained us- 
ing the theory presented above for two core-halo systems 
characterized by different values of initial r m and i> m , and 
compare these distributions with the ones obtained using 
the dynamics simulations. Excellent agreement is found 
in all cases. In Fig. [3] we also present the distribution 
functions obtained using the violent relaxation theory, 
Eq. (|4]). It is clear that this theory is unable to describe 
relaxation of oscillating plasmas. 

Up to now we have considered plasmas which at t = 
where uniformly distributed. This, however, is not very 
usual in practice and more realistically one might expect 
a initially thermallized distribution of the form 

/o(r,v) = - 2 l 22 e(r m -r)e-& . (11) 

The procedure is then to discretize the Gaussian part of 
this distribution into p levels. At the lowest order, we 
can take p — 1 and approximate Eq. (fTTj) by Eq. ([8]). 
To have equal energy, both distributions must have the 
same values of (v 2 ). This requires that v m = la. The 



final relaxed distribution of this core-halo system should 
then be given by Eq. (fTOf with r\\ = l/47r 2 r^cr 2 . In 
Fig. 4 we plot N(r)/N}, and N(v)/Nb and compare these 
distributions with the ones obtained using the dynamics 
simulation in which the initial particle positions and ve- 
locities were distributed according to Eq. (jTTJ) . In spite of 
the crudeness of the one level approximation, an amaz- 
ingly good agreement between the theory and the sim- 
ulations is obtained without any adjustable parameters. 
We have checked that this good agreement persists for 
other values of a and r m . 

In this Letter we have studied confined one component 
plasmas of charges interacting by unscreened Coulomb 
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FIG. 4: (a) Position and (b) velocity distributions. Solid 
curves are the theoretical predictions and points are the result 
of dynamics simulations. Initial t = distribution is a thermal 
one given by Eq. JTT]) r m = and a = 0.64^. 



potential. Unlike normal gases with short range forces, 
non-neutral plasmas do not evolve to the state of ther- 
modynamic equilibrium. Instead collisionless relaxation 
culminates in a stationary state in which the detailed bal- 
ance is violated. Using a combination of non-equilibrium 
statistical mechanics and the theory of parametric reso- 
nances it is, nevertheless, possible to a priory predict the 
distribution functions for the final stationary state. Un- 
like the normal thermodynamic equilibrium, this state, 
however, explicitly depends on the initial distribution of 
particle velocities and positions. 
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